Why mapping in R?
Spatial data & Coordinate reference system
Importation and attribute manipulation of vector data with
sfImportation of raster data manipulation with
rasterGeometry manipulations for vector and raster data
Examples of maps created with R
Questions, discussion, and use of your data
30 avril 2018
Workshop outline
Why?
Mapping in Ecology?
show where your plots are
show how variables are distributed spatially
- show results of spatial analyses
Why using R as a GIS?
- Open-source, free
- Benefit from a very active community
- Very large number of packages
Why using R as a GIS?
Open-source, free
- Workflow and replicability
- import, format, analyze, visualize, export your data
- repeat with new data
- create your own functions/packages
Why using R as a GIS?
Open-source, free
Workflow and replicability
- Quite efficient
- well-defined spatial classes
- can read/write/convert many formats
Why using R as a GIS?
Open-source, free
Workflow and replicability
Quite efficient
- Many interfaces to other software
But…
Possible limitations - Geo-referencing - Visualizing large spatial objects - Watershed analysis - …
Spatial data
Vector data
Raster data
Coordinate reference system
Coordinate reference system (CRS)
- A place on the earth is specified by a latitude and longitude or X/Y coordinates
- The coordinates are based on a mathematical model of the shape of the earth
- Some mathematical formulas involve a transformation from the 3-D globle to a 2-D map
Projected vs geographic CRS
An geographic (or unprojected) CRS uses latitude and longitude coordinates and references the earth as a 3-D object
A projected CRS uses X and Y coordinates as a 2-D representation of the earth
Define CRS with EPSG or proj4string
Your geographic files have a CRS - but it’s not always defined - Vector and raster spatial data was created based on a specific CRS - Usually spatial files have metadata that tells you the CRS
You can use either an EPSG or proj4string to define the CRS.
- The EPSG code is a numeric representation of a CRS (e.g, 4326)
- The proj4string is a full set of parameters spelled out in a string (e.g., “+proj=longlat +ellps=WGS84 +no_defs”)
To find CRS in any format Spatial Reference
Vector data
Intro to Simple Features
Why use the sf package when sp is already tried and tested?
Simple features refers to a formal standard (ISO 19125-1:2004) that describes how objects in the real world can be represented in computers
Soon to be the successor of
sp(spdeveloped in 2003,sfintroduced in 2016)sfincorporates the functionality of the 3 main packages of thespparadigm in a single, cohesive whole:spfor the class system;rgdalfor reading and writing data;rgeosfor spatial operations undertaken by GEOS).
Intro to Simple Features
Why use the sf package when sp is already tried and tested?
sfobjects are easy to manipulate- Spatial objects are stored as data frames, with the feature geometries stored in list-columns
- All functions begin with
st_for easy tab completion - Functions are pipe-friendly
%>% dplyrandtidyrverbs have been defined for the sf objects- Fast reading and writing of data
ggplotfriendly
- GREAT documentation! sf vignette
Intro to Simple Features
Import your dataset in R
- Non-spatial data
- Water quality measures in Montreal 1
- Vector data
- Points - sample points of water quality measures in Montreal 1
- Lines - streams and rivers of Montreal 1
- Polygons - polygons of land use types 2
- Raster data
- Raster - canopy index of Montreal 2
- Raster - altitude (retrieved directly from R)
Packages
library(sf) # for simple features vector library(lwgeom) # for st_make_valid library(dplyr) # for data manipulation library(RColorBrewer) # for color palette library(raster) # for raster data library(mapview) # for interactive map
MULTIPOINT
Load and import sample points from a csv
# Create a new directory to download data
if(!dir.exists("data")) dir.create("data")
# Download csv file from web page in your working directory
if (!file.exists("data/ruisso.csv"))
download.file("http://donnees.ville.montreal.qc.ca/dataset/86843d31-4251-4002-b10d-620aaa751092/resource/adad6c48-fb48-40fc-a031-1ac870d978b4/download/scri03.-infor03.02-informatique03.02.07-donneesouvertesrsmaruissostationsstations_ruisso.csv",
destfile = "data/ruisso.csv")
Load and import sample points from a csv
# Read csv file in R
ruisso <- read.csv("data/ruisso.csv", header = T, dec = ",")
str(ruisso)
## 'data.frame': 66 obs. of 7 variables: ## $ Plan.d.eau : Factor w/ 31 levels "Bassin de La Brunante",..: 20 20 20 20 20 20 20 20 20 20 ... ## $ Point.d.échantillonnage: Factor w/ 66 levels "AAO-0.0","AAO-1.5P1",..: 1 2 3 4 5 6 7 8 9 10 ... ## $ Localisation : Factor w/ 65 levels "200m à l'est du Boul.Morgan et 280m au sud de la transcanadienne, affluent provenant de Ste-Anne-de-Bellevue et de la partie O "| __truncated__,..: 49 51 64 34 63 13 12 1 10 11 ... ## $ Administration : Factor w/ 21 levels "Ahuntsic-Cartierville",..: 14 14 18 7 18 3 4 3 3 3 ... ## $ LATITUDE : num 45.5 45.4 45.4 45.4 45.4 ... ## $ LONGITUDE : num -73.9 -73.9 -73.9 -73.9 -73.9 ... ## $ Activité : Factor w/ 2 levels "Actif","Inactif": 1 2 2 1 1 1 2 2 1 1 ...
Convert to sf MULTIPOINT object
ruisso_sf <- st_as_sf(
x = ruisso,
coords = c("LONGITUDE", "LATITUDE"),
crs = 4326)
ruisso_sf
## Simple feature collection with 66 features and 5 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -73.93704 ymin: 45.42516 xmax: -73.50467 ymax: 45.69734 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## First 10 features: ## Plan.d.eau Point.d.échantillonnage ## 1 Rivière à l'Orme AAO-0.0 ## 2 Rivière à l'Orme AAO-1.5P1 ## 3 Rivière à l'Orme AAO-2.0P4 ## 4 Rivière à l'Orme AAO-3.3P6 ## 5 Rivière à l'Orme AAO-3.5 ## 6 Rivière à l'Orme AAO-3.6 ## 7 Rivière à l'Orme AAO-3.9P7 ## 8 Rivière à l'Orme AAO-6.4 ## 9 Rivière à l'Orme AAO-6.4P12 ## 10 Rivière à l'Orme AAO-6.5 ## Localisation ## 1 Pierrefonds, boul. Gouin O, 40m au nord de la rue de l'Anse à l'Orme, exutoire au lac des Deux Montagnes. ## 2 Pierrefonds, N ponceau boul.Gouin, 1500m en amont exutoire, branche provenant de l'est. ## 3 Ste-A.-de-Bellevue, branche drainant secteur ouest, 140m à l'est de la rue Leslie Dowker. ## 4 Kirkland, 60m au sud de l'intersection des rues de l'Anse à l'Orme et de Timberley trail, derrière le dépôt à neige. ## 5 Sainte-Anne-de-Bellevue, 10m au nord du ch. Ste-Marie, 200m à l'ouest du ch. Anse à l'Orme. ## 6 Beaconsfield, 250m à l'est de la rue Lee et 25m au sud de l'autoroute 40, en amont du pluvial. ## 7 Beaconsfield, 240m à l'est de la rue Lee et 400m au sud de l'A40, embranchement provenant de la zone boisée entourant le boul. Lakeview. ## 8 200m à l'est du Boul.Morgan et 280m au sud de la transcanadienne, affluent provenant de Ste-Anne-de-Bellevue et de la partie O et S de la zone industriel de Baie d'Urfée. ## 9 Baie d'Urfée, boul. Morgan côté est, 250m au sud de l'autoroute 40, affluent provenant des zones résidentilelles de Baie d'Urfé. ## 10 Baie d'Urfée, boul.Morgan côté ouest, 250m au sud de l'A40. ## Administration Activité geometry ## 1 Pierrefonds-Roxboro Actif POINT (-73.93704 45.45022) ## 2 Pierrefonds-Roxboro Inactif POINT (-73.91931 45.44744) ## 3 Sainte-Anne-de-Bellevue Inactif POINT (-73.91535 45.44288) ## 4 Kirkland Actif POINT (-73.90147 45.43689) ## 5 Sainte-Anne-de-Bellevue Actif POINT (-73.90144 45.43566) ## 6 Baie d'Urfé Actif POINT (-73.9012 45.43462) ## 7 Beaconsfield Inactif POINT (-73.9002 45.43105) ## 8 Baie d'Urfé Inactif POINT (-73.91347 45.42674) ## 9 Baie d'Urfé Actif POINT (-73.91549 45.42531) ## 10 Baie d'Urfé Actif POINT (-73.91565 45.42542)
Simple mapping of MULTIPOINT
Instead of creating a single map, as with sp object, the default plot of sf object creates multiple maps, one for each attribute, which can be useful for exploring the spatial distribution of different variables.
plot(ruisso_sf)
Simple mapping of MULTIPOINT
To plot only the geometry and not all attributes, we retrieve the geometry list-column using st_geometry():
plot(st_geometry(ruisso_sf))
Load and import sample data from a csv
# Download csv file from web page in your working directory
if (!file.exists("data/donnees_ruisso_2016.csv"))
download.file("http://donnees.ville.montreal.qc.ca/dataset/8c149ace-7b2f-4041-99ec-3bdbef5dcee6/resource/38c8eb26-46a1-4fc2-87a0-5c88e94cee8e/download/donnees_ruisso_2016.csv",
destfile = "data/ruisso_data.csv")
Load and import sample data from a csv
# Read csv file in R
ruisso_data <- read.csv("data/ruisso_data.csv", header = T, dec = ",")
str(ruisso_data)
## 'data.frame': 345 obs. of 40 variables: ## $ Point.d.échantillonnage: Factor w/ 50 levels "AAO-0.0","AAO-3.3P6",..: 1 1 1 1 1 1 1 2 2 2 ... ## $ Date : Factor w/ 32 levels "2016/05/09","2016/05/12",..: 2 6 9 13 18 23 30 2 6 9 ... ## $ Raison.d.annulation : Factor w/ 4 levels "","À sec","Niveau trop faible",..: 1 1 1 1 1 1 1 1 1 1 ... ## $ X.OD : num 110 74 72 79 79 76.8 96.9 110 81 93 ... ## $ O2..mg.L. : num 10.7 7.4 6.7 6.8 7.3 7.1 11 11.6 8.1 9.2 ... ## $ Conductivité..µs.cm2. : int 514 1434 1474 1346 793 1286 1414 1589 1639 897 ... ## $ pH : num 7.9 7.9 7.7 7.9 7.8 7.7 7.7 7.9 8 7.9 ... ## $ Température..oC. : num 16.9 15 19.1 22.3 19.1 18.8 9.4 12.9 15.6 15.6 ... ## $ Signe.COLI : Factor w/ 4 levels "","<","=",">": 2 3 3 3 3 3 3 3 3 3 ... ## $ COLI : int 10 54 230 300 550 72 99 1200 4400 33000 ... ## $ MÉTÉO : int 1 1 0 1 -1 -2 -2 1 1 0 ... ## $ Ag..µg.L. : num 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ... ## $ Al..µg.L. : num 214 60 761 214 353 207 259 81 95 2600 ... ## $ As..µg.L. : num 0.4 0.5 0.9 0.8 0.5 0.5 0.4 0.3 0.4 1.1 ... ## $ Ba..µg.L. : num 37 69 58 67 56 61 77 60 56 63 ... ## $ Be..µg.L. : num 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ... ## $ Ca..µg.L. : int 44200 109000 104000 98200 71200 102000 122000 125000 123000 72700 ... ## $ Cd..µg.L. : num 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ... ## $ Co..µg.L. : num 0.2 0.2 0.8 0.3 0.4 0.2 0.3 0.2 0.2 2.7 ... ## $ COT..µg.L. : num 7.2 6.4 6.1 15.1 4.3 5.3 4.2 3.8 3.7 23.6 ... ## $ Cr..µg.L. : num 0.6 0.3 2.2 0.5 1.3 0.7 0.9 0.3 0.3 9.9 ... ## $ Cu..µg.L. : num 1.9 3.2 4.3 1.6 3 3.2 3.1 3.3 8 30 ... ## $ Fe..µg.L. : int 364 271 1200 393 533 367 422 237 258 3950 ... ## $ K..µg.L. : int 1970 4190 4600 4180 3150 4100 5590 4070 4840 4790 ... ## $ Mg..µg.L. : int 15600 37100 39300 36000 19500 33700 35900 36300 35400 20900 ... ## $ Mn..µg.L. : num 40.4 70.3 82.8 65.1 44.2 18.9 24.7 57.2 68.3 261 ... ## $ Mo..µg.L. : num 1.5 4.1 4.1 3.9 2.9 3.3 3.8 2.3 3.1 2.6 ... ## $ Na..µg.L. : int 46600 100000 100000 100000 58200 100000 100000 100000 100000 92300 ... ## $ NH3..µg.L. : num 30 160 290 90 120 40 70 60 100 490 ... ## $ Ni..µg.L. : num 1.2 2.3 3 2 2.3 2.7 2.4 2 2.1 7.6 ... ## $ Ptot..µg.L. : num 37 33 129 67 55 31 36 23 207 430 ... ## $ Pb..µg.L. : num 0.5 0.2 3.2 0.9 0.8 0.7 0.6 0.2 0.2 5.6 ... ## $ MES...mg.L. : num 6.8 5.7 35.3 9 10.3 3.6 5.8 3.4 5 161 ... ## $ Sb..µg.L. : num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1.2 ... ## $ Se..µg.L. : num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ... ## $ Sn2.ug.L : num 1 1 1 1 1 1 1 1 1 1 ... ## $ Tl.ug.L : num 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ... ## $ U..µg.L. : num 0.8 1.7 1.8 1.8 1.2 1.3 2.2 1.9 1.6 0.9 ... ## $ V..µg.L. : num 0.9 0.6 2.8 1.7 1.7 1.2 1.2 0.7 1 9.4 ... ## $ Zn..µg.L. : num 7 7 12 7 7 7 7 7 7 108 ...
Attribute manipulations
Useful links for sf manipulations
Join these sampled data as attributes
Combining data from different sources is a common task in data preparation. left_join() from dplyr do this by combining tables based on a shared ‘key’ variable.
See dplyr and tidyr cheatsheet for other join functions.
ruisso_sf <- left_join(ruisso_sf, ruisso_data, by = "Point.d.échantillonnage") names(ruisso_sf)
## [1] "Plan.d.eau" "Point.d.échantillonnage" ## [3] "Localisation" "Administration" ## [5] "Activité" "Date" ## [7] "Raison.d.annulation" "X.OD" ## [9] "O2..mg.L." "Conductivité..µs.cm2." ## [11] "pH" "Température..oC." ## [13] "Signe.COLI" "COLI" ## [15] "MÉTÉO" "Ag..µg.L." ## [17] "Al..µg.L." "As..µg.L." ## [19] "Ba..µg.L." "Be..µg.L." ## [21] "Ca..µg.L." "Cd..µg.L." ## [23] "Co..µg.L." "COT..µg.L." ## [25] "Cr..µg.L." "Cu..µg.L." ## [27] "Fe..µg.L." "K..µg.L." ## [29] "Mg..µg.L." "Mn..µg.L." ## [31] "Mo..µg.L." "Na..µg.L." ## [33] "NH3..µg.L." "Ni..µg.L." ## [35] "Ptot..µg.L." "Pb..µg.L." ## [37] "MES...mg.L." "Sb..µg.L." ## [39] "Se..µg.L." "Sn2.ug.L" ## [41] "Tl.ug.L" "U..µg.L." ## [43] "V..µg.L." "Zn..µg.L." ## [45] "geometry"
Data cleaning
Keep only active sites
ruisso_sf1 <- filter(ruisso_sf, Activité == "Actif") # same as # ruisso_sf1 <- filter(ruisso_sf, Activité != "Inactif") dim(ruisso_sf)
## [1] 361 45
dim(ruisso_sf1)
## [1] 345 45
Data cleaning
Remove unwanted variables
ruisso_sf2 <- dplyr::select(ruisso_sf1,
-Localisation,
-Activité,
-Raison.d.annulation,
-Signe.COLI,
-MÉTÉO)
Note: Both raster and dplyr packages have a function called select(). To avoid an error message when both packages are loaded, we use the long-form function name: dplyr::select().
Data cleaning
names(ruisso_sf2)
## [1] "Plan.d.eau" "Point.d.échantillonnage" ## [3] "Administration" "Date" ## [5] "X.OD" "O2..mg.L." ## [7] "Conductivité..µs.cm2." "pH" ## [9] "Température..oC." "COLI" ## [11] "Ag..µg.L." "Al..µg.L." ## [13] "As..µg.L." "Ba..µg.L." ## [15] "Be..µg.L." "Ca..µg.L." ## [17] "Cd..µg.L." "Co..µg.L." ## [19] "COT..µg.L." "Cr..µg.L." ## [21] "Cu..µg.L." "Fe..µg.L." ## [23] "K..µg.L." "Mg..µg.L." ## [25] "Mn..µg.L." "Mo..µg.L." ## [27] "Na..µg.L." "NH3..µg.L." ## [29] "Ni..µg.L." "Ptot..µg.L." ## [31] "Pb..µg.L." "MES...mg.L." ## [33] "Sb..µg.L." "Se..µg.L." ## [35] "Sn2.ug.L" "Tl.ug.L" ## [37] "U..µg.L." "V..µg.L." ## [39] "Zn..µg.L." "geometry"
Data cleaning
Rename some variables
ruisso_sf3 <- rename(ruisso_sf2, river = Plan.d.eau, sample_pts = Point.d.échantillonnage, dissolve_O = X.OD) names(ruisso_sf3)
## [1] "river" "sample_pts" ## [3] "Administration" "Date" ## [5] "dissolve_O" "O2..mg.L." ## [7] "Conductivité..µs.cm2." "pH" ## [9] "Température..oC." "COLI" ## [11] "Ag..µg.L." "Al..µg.L." ## [13] "As..µg.L." "Ba..µg.L." ## [15] "Be..µg.L." "Ca..µg.L." ## [17] "Cd..µg.L." "Co..µg.L." ## [19] "COT..µg.L." "Cr..µg.L." ## [21] "Cu..µg.L." "Fe..µg.L." ## [23] "K..µg.L." "Mg..µg.L." ## [25] "Mn..µg.L." "Mo..µg.L." ## [27] "Na..µg.L." "NH3..µg.L." ## [29] "Ni..µg.L." "Ptot..µg.L." ## [31] "Pb..µg.L." "MES...mg.L." ## [33] "Sb..µg.L." "Se..µg.L." ## [35] "Sn2.ug.L" "Tl.ug.L" ## [37] "U..µg.L." "V..µg.L." ## [39] "Zn..µg.L." "geometry"
Data cleaning
Remove symbole from column names
names(ruisso_sf3) <- gsub("\\..*", "", names(ruisso_sf3))
names(ruisso_sf3)
## [1] "river" "sample_pts" "Administration" "Date" ## [5] "dissolve_O" "O2" "Conductivité" "pH" ## [9] "Température" "COLI" "Ag" "Al" ## [13] "As" "Ba" "Be" "Ca" ## [17] "Cd" "Co" "COT" "Cr" ## [21] "Cu" "Fe" "K" "Mg" ## [25] "Mn" "Mo" "Na" "NH3" ## [29] "Ni" "Ptot" "Pb" "MES" ## [33] "Sb" "Se" "Sn2" "Tl" ## [37] "U" "V" "Zn" "geometry"
Data cleaning
Remove accent from column names
names(ruisso_sf3) <- gsub("é", "e", names(ruisso_sf3))
names(ruisso_sf3)
## [1] "river" "sample_pts" "Administration" "Date" ## [5] "dissolve_O" "O2" "Conductivite" "pH" ## [9] "Temperature" "COLI" "Ag" "Al" ## [13] "As" "Ba" "Be" "Ca" ## [17] "Cd" "Co" "COT" "Cr" ## [21] "Cu" "Fe" "K" "Mg" ## [25] "Mn" "Mo" "Na" "NH3" ## [29] "Ni" "Ptot" "Pb" "MES" ## [33] "Sb" "Se" "Sn2" "Tl" ## [37] "U" "V" "Zn" "geometry"
Summarise attributes by sample plots
There are multiple measurements during the summer. We will use the mean of water quality measurements.
ruisso_sf3
## Simple feature collection with 345 features and 39 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -73.93704 ymin: 45.42516 xmax: -73.50467 ymax: 45.69734 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## First 10 features: ## river sample_pts Administration Date dissolve_O ## 1 Rivière à l'Orme AAO-0.0 Pierrefonds-Roxboro 2016/05/12 110.0 ## 2 Rivière à l'Orme AAO-0.0 Pierrefonds-Roxboro 2016/06/01 74.0 ## 3 Rivière à l'Orme AAO-0.0 Pierrefonds-Roxboro 2016/06/22 72.0 ## 4 Rivière à l'Orme AAO-0.0 Pierrefonds-Roxboro 2016/08/02 79.0 ## 5 Rivière à l'Orme AAO-0.0 Pierrefonds-Roxboro 2016/08/22 79.0 ## 6 Rivière à l'Orme AAO-0.0 Pierrefonds-Roxboro 2016/09/20 76.8 ## 7 Rivière à l'Orme AAO-0.0 Pierrefonds-Roxboro 2016/10/25 96.9 ## 8 Rivière à l'Orme AAO-3.3P6 Kirkland 2016/05/12 110.0 ## 9 Rivière à l'Orme AAO-3.3P6 Kirkland 2016/06/01 81.0 ## 10 Rivière à l'Orme AAO-3.3P6 Kirkland 2016/06/22 93.0 ## O2 Conductivite pH Temperature COLI Ag Al As Ba Be Ca Cd ## 1 10.7 514 7.9 16.9 10 0.1 214 0.4 37 0.1 44200 0.1 ## 2 7.4 1434 7.9 15.0 54 0.1 60 0.5 69 0.1 109000 0.1 ## 3 6.7 1474 7.7 19.1 230 0.1 761 0.9 58 0.1 104000 0.1 ## 4 6.8 1346 7.9 22.3 300 0.1 214 0.8 67 0.1 98200 0.1 ## 5 7.3 793 7.8 19.1 550 0.1 353 0.5 56 0.1 71200 0.1 ## 6 7.1 1286 7.7 18.8 72 0.1 207 0.5 61 0.1 102000 0.1 ## 7 11.0 1414 7.7 9.4 99 0.1 259 0.4 77 0.1 122000 0.1 ## 8 11.6 1589 7.9 12.9 1200 0.1 81 0.3 60 0.1 125000 0.1 ## 9 8.1 1639 8.0 15.6 4400 0.1 95 0.4 56 0.1 123000 0.1 ## 10 9.2 897 7.9 15.6 33000 0.1 2600 1.1 63 0.1 72700 0.1 ## Co COT Cr Cu Fe K Mg Mn Mo Na NH3 Ni Ptot Pb ## 1 0.2 7.2 0.6 1.9 364 1970 15600 40.4 1.5 46600 30 1.2 37 0.5 ## 2 0.2 6.4 0.3 3.2 271 4190 37100 70.3 4.1 100000 160 2.3 33 0.2 ## 3 0.8 6.1 2.2 4.3 1200 4600 39300 82.8 4.1 100000 290 3.0 129 3.2 ## 4 0.3 15.1 0.5 1.6 393 4180 36000 65.1 3.9 100000 90 2.0 67 0.9 ## 5 0.4 4.3 1.3 3.0 533 3150 19500 44.2 2.9 58200 120 2.3 55 0.8 ## 6 0.2 5.3 0.7 3.2 367 4100 33700 18.9 3.3 100000 40 2.7 31 0.7 ## 7 0.3 4.2 0.9 3.1 422 5590 35900 24.7 3.8 100000 70 2.4 36 0.6 ## 8 0.2 3.8 0.3 3.3 237 4070 36300 57.2 2.3 100000 60 2.0 23 0.2 ## 9 0.2 3.7 0.3 8.0 258 4840 35400 68.3 3.1 100000 100 2.1 207 0.2 ## 10 2.7 23.6 9.9 30.0 3950 4790 20900 261.0 2.6 92300 490 7.6 430 5.6 ## MES Sb Se Sn2 Tl U V Zn geometry ## 1 6.8 0.5 0.5 1 0.2 0.8 0.9 7 POINT (-73.93704 45.45022) ## 2 5.7 0.5 0.5 1 0.2 1.7 0.6 7 POINT (-73.93704 45.45022) ## 3 35.3 0.5 0.5 1 0.2 1.8 2.8 12 POINT (-73.93704 45.45022) ## 4 9.0 0.5 0.5 1 0.2 1.8 1.7 7 POINT (-73.93704 45.45022) ## 5 10.3 0.5 0.5 1 0.2 1.2 1.7 7 POINT (-73.93704 45.45022) ## 6 3.6 0.5 0.5 1 0.2 1.3 1.2 7 POINT (-73.93704 45.45022) ## 7 5.8 0.5 0.5 1 0.2 2.2 1.2 7 POINT (-73.93704 45.45022) ## 8 3.4 0.5 0.5 1 0.2 1.9 0.7 7 POINT (-73.90147 45.43689) ## 9 5.0 0.5 0.5 1 0.2 1.6 1.0 7 POINT (-73.90147 45.43689) ## 10 161.0 1.2 0.5 1 0.2 0.9 9.4 108 POINT (-73.90147 45.43689)
Summarise attributes by sample plots
We will use the mean of water quality measurements.
ruisso_mean <- ruisso_sf3 %>% group_by(sample_pts) %>% summarise_at(vars(O2:Zn), mean, na.rm = TRUE) ruisso_mean
## Simple feature collection with 50 features and 35 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -73.93704 ymin: 45.42516 xmax: -73.50467 ymax: 45.69734 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## # A tibble: 50 x 36 ## sample_pts O2 Conductivite pH Temperature COLI Ag Al ## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 AAO-0.0 8.14 1180. 7.80 17.2 188. 0.100 295. ## 2 AAO-3.3P6 9.50 1378. 7.89 16.4 18599. 0.100 432. ## 3 AAO-3.5 8.79 1281. 7.74 14.3 612. 0.100 150. ## 4 AAO-3.6 9.00 1128. 7.97 16.1 461. 0.100 217. ## 5 AAO-6.4P12 10.1 1461. 7.96 18.2 147. 0.100 33.1 ## 6 AAO-6.5 7.40 567. 8.12 16.5 1219. 0.100 109. ## 7 ADM-1 9.97 333. 8.19 21.2 58.6 0.100 119. ## 8 ANG-2 9.70 399. 8.33 20.2 41.0 0.100 60.7 ## 9 BER-0.0 7.18 803. 7.61 17.1 305. 0.100 140. ## 10 BER-0.7P1 9.25 751. 7.76 17.1 1431. 0.100 75.7 ## # ... with 40 more rows, and 28 more variables: As <dbl>, Ba <dbl>, ## # Be <dbl>, Ca <dbl>, Cd <dbl>, Co <dbl>, COT <dbl>, Cr <dbl>, Cu <dbl>, ## # Fe <dbl>, K <dbl>, Mg <dbl>, Mn <dbl>, Mo <dbl>, Na <dbl>, NH3 <dbl>, ## # Ni <dbl>, Ptot <dbl>, Pb <dbl>, MES <dbl>, Sb <dbl>, Se <dbl>, ## # Sn2 <dbl>, Tl <dbl>, U <dbl>, V <dbl>, Zn <dbl>, geometry <POINT [°]>
Simple mapping of attributes
plot(ruisso_mean)
Simple mapping of attributes
myPal <- colorRampPalette(c("blue", "red"))
plot(ruisso_mean["Temperature"],
pal = myPal, nbreaks = 30, pch = 19, key.pos = 1,
main = "Water temperature in streams of Montreal")
Export your MULTIPOINT to a Shapefile
We can write a simple features object to a file (e.g. a shapefile) using the st_write() function in sf, which needs at least two arguments, the object and a filename:
st_write(ruisso_mean, "data/ruisso.shp", delete_dsn = T)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ## ESRI Shapefile driver
Load and import land use polygons
We will now import a land use vector map of Montreal.
Load and import land use polygons
Because the original dataset was composed of multiple individual shapefiles and were long to download, some manipulations were performed before (download from webpage, combine shapefiles together, polygon validation). You can read a ready-to-use shapefile directly from
Don’t run these lines now!
# Download shapefiles
download.file("http://cmm.qc.ca/fileadmin/user_upload/geomatique/UtilisationDuSol/2016_Shapefiles/660-US-2016.zip", destfile = "data/landuse.zip")
# Unzip the main folder and name it "landuse"
unzip("landuse.zip", exdir="data/landuse")
# get all the zip files inside the main folder "landuse"
zipF <- list.files(path = "data/landuse/", pattern = "*.zip", full.names = TRUE)
# unzip all your files
plyr::ldply(.data = zipF, .fun = unzip, exdir = "data/landuse")
# Get the names of the land use shapefiles from the folder "landuse"
shp_files <- list.files(path = "data/landuse", pattern = ".shp")
shp_files <- sub(".shp", "", shp_files)
Load and import land use polygons
Don’t run these lines now!
# Read all the shapefiles
LU <- list()
for(f in shp_files) {
LU[[f]] <- st_read(dsn = "data/landuse/", layer = f)
}
# Combine all shapefiles together
LU_all <- do.call(rbind, LU)
# Write as a GeoPackage
st_write(LU_all, "data/LU_all.gpkg", driver = "GPKG")
Load and import land use polygons
# Read GeoPackage in R LU <- st_read(dsn = "data/LU_all.gpkg")
## Reading layer `LU_all' from data source `/Users/mariehbrice/Documents/GitHub/Rspatial/data/LU_all.gpkg' using driver `GPKG' ## Simple feature collection with 107233 features and 38 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: 265961.2 ymin: 5027324 xmax: 306814.6 ymax: 5063078 ## epsg (SRID): NA ## proj4string: +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Notice that the EPSG code is not defined
Defining the crs with st_set_crs()
To define the CRS when it’s missing we can use st_set_crs(). However, replacing crs does not reproject data. If we want to transform coordinate and reproject them, we would use st_transform().
LU <- st_set_crs(LU, 32188)
Mapping MULTIPOLYGON
We can map only one land use at a time by subsetting first the sf object. UTIL_SOL==1000 represents water.
plot(st_geometry(subset(LU, UTIL_SOL==1000)))
Corrupt or invalid geometries?
Geometry validity refers to essential properties of polygons, such as non-self intersecting, holes being inside polygons, less than 4 points, non-closing segments. Invalid geometry can be ignored for mapping but cause problem for spatial operations.
st_is_valid() returns a logical vector indicating for each polygon geometries whether it is topologically valid:
(invalid_poly <- which(!st_is_valid(LU)))
## [1] 862 1149 1251 1406 2071 2119 3479 3837 4838 5586 ## [11] 7908 7975 8057 8073 8111 8132 8501 8980 9001 9149 ## [21] 9252 9638 9715 10381 10424 10910 10929 11342 11801 12094 ## [31] 13135 14420 16639 20129 20228 20678 21322 21445 22679 22800 ## [41] 23330 23684 26948 26951 27117 27122 27144 27152 27479 28796 ## [51] 28917 28975 29932 30268 31097 31945 32140 32475 32543 32557 ## [61] 32640 32946 33281 33334 33834 35010 36204 36865 38359 39235 ## [71] 39251 39255 39276 39535 39549 39919 41952 42034 42085 42439 ## [81] 42484 42781 43534 44540 45171 50868 56315 58126 58127 59198 ## [91] 59219 64288 66720 68063 68539 68639 68701 68744 69507 69977 ## [101] 71081 78557 78740 79890 81941 83127 87246 87384 87419 87434 ## [111] 93149 93722 93730 93899 93980 94074 94159 94593 94839 99756 ## [121] 100103 101174 101228 101260 101330 103278 103877 104673 105625 105951 ## [131] 106710 106874 106900 106945 107213
Using the argument reason = TRUE returns the reason for invalidity:
st_is_valid(LU[invalid_poly,], reason = TRUE)[1:8]
## [1] "Nested shells[295499.0604 5045532.3412]" ## [2] "Nested shells[298395.5742 5045168.6803]" ## [3] "Self-intersection[297511.205 5032954.7058]" ## [4] "Nested shells[273102.34356 5035592.911633]" ## [5] "Nested shells[274053.456247 5036041.908298]" ## [6] "Nested shells[288581.300144 5036823.875472]" ## [7] "Self-intersection[294414.9683 5031596.2768]" ## [8] "Self-intersection[295265.9577 5045076.9297]"
Corrupt or invalid geometries?
par(mfrow = c(1,2)) plot(st_geometry(LU[862,]), main = "Nested shells") plot(st_geometry(LU[1251,]), main = "Self-intersection")
Corrupt or invalid geometries?
We can use st_make_valid() from lwgeom to make an invalid geometry valid
LU_val <- st_make_valid(LU) which(!st_is_valid(LU_val)) # yeah!
## integer(0)
Exercise
Exercise
- Import a shapefile of watercourses in Montreal (courseau.shp) using
st_read(dsn = path_to_file, layer = file_name, driver = "ESRI Shapefile")and name itcourseauIf you did not load the shapefile before the workshop you can download it from this link using
download.file(): http://donnees.ville.montreal.qc.ca/dataset/c128aff5-325c-4599-ab66-1c9d0b3abc94/resource/a37e11d4-f0a3-46a7-8636-76754fad72b3/download/courseau.zip Create a simple map of the geometry
Create a simple thematic map of watercourse
TYPE. You can change the default colors using the argumentpal.
## Reading layer `courseau' from data source `/Users/mariehbrice/Documents/GitHub/Rspatial/data/courseau' using driver `ESRI Shapefile' ## Simple feature collection with 1306 features and 5 fields ## geometry type: MULTILINESTRING ## dimension: XY ## bbox: xmin: -73.97268 ymin: 45.41593 xmax: -73.4975 ymax: 45.69939 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs
Exercise
## Simple feature collection with 1306 features and 5 fields ## geometry type: MULTILINESTRING ## dimension: XY ## bbox: xmin: -73.97268 ymin: 45.41593 xmax: -73.4975 ymax: 45.69939 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## First 10 features: ## OBJECTID_1 NOM TYPE Shape_Le_1 NuméroRui ## 1 1 rivière à l'Orme rivière 177.95299 <NA> ## 2 2 rivière à l'Orme rivière 128.51146 <NA> ## 3 3 <NA> fossé 172.42988 <NA> ## 4 4 rivière à l'Orme rivière 216.66838 <NA> ## 5 8 <NA> fossé 540.29539 <NA> ## 6 9 <NA> ruisseau 97.66412 <NA> ## 7 10 <NA> ruisseau 162.30584 <NA> ## 8 11 <NA> ruisseau 210.75794 <NA> ## 9 14 <NA> ruisseau 608.28651 <NA> ## 10 16 <NA> fossé 267.61108 <NA> ## geometry ## 1 MULTILINESTRING ((-73.9107 ... ## 2 MULTILINESTRING ((-73.90824... ## 3 MULTILINESTRING ((-73.90667... ## 4 MULTILINESTRING ((-73.91472... ## 5 MULTILINESTRING ((-73.91029... ## 6 MULTILINESTRING ((-73.9206 ... ## 7 MULTILINESTRING ((-73.91969... ## 8 MULTILINESTRING ((-73.91727... ## 9 MULTILINESTRING ((-73.91727... ## 10 MULTILINESTRING ((-73.91212...
Exercise
Exercise
Raster data
Useful links
raster classes
The R package raster provides three main classes of raster object (more details here):
RasterLayer imports a single-layer (variable) raster
RasterStack imports in one single object several single-layer (variable) rasters stored in one or different files
RasterBrick imports in one single object several single-layer (variable) rasters stored in one single file.
Import raster data
We now import raster data use a .tif file of a canopy index from Montreal.
# Download tif file from web page in your working directory
if (!file.exists("data/canopee.zip"))
download.file("http://cmm.qc.ca/fileadmin/user_upload/geomatique/IndiceCanopee/2015/660_Canopee2015_3m.zip",
destfile = "data/canopee.zip")
unzip("data/canopee.zip", exdir = "data")
# Read tif in R using raster
# The file named "660_CLASS_3m.tif" contains the canopy index for all the Montreal area, so we can read this file only
canopee_mtl <- raster("data/660_CLASS_3m.tif")
raster
The canopy index raster has values from 1 to 5, has 35754 pixels by row and 40854 pixels by column and resolution of 1m x 1m.
canopee_mtl
## class : RasterLayer ## dimensions : 35754, 40854, 1460693916 (nrow, ncol, ncell) ## resolution : 1, 1 (x, y) ## extent : 265961, 306815, 5027324, 5063078 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ## data source : /Users/mariehbrice/Documents/GitHub/Rspatial/data/660_CLASS_3m.tif ## names : X660_CLASS_3m ## values : 1, 5 (min, max)
Simple mapping of raster
Similar to the sf package, raster also provides plot methods for its own classes.
plot(canopee_mtl)
Retrieving free GIS data: getData
In package raster, getData() function generates requests to access to different spatial datasets (altitude, climate, country boundary). Argument name specifies the dataset you wish to download.
head(getData("ISO3"))
## ISO3 NAME ## 1 ABW Aruba ## 2 AFG Afghanistan ## 3 AGO Angola ## 4 AIA Anguilla ## 5 ALA Åland ## 6 ALB Albania
Retrieving free GIS data: getData
Retrieve a raster of altitude for Canada.
altCAN <- getData(name = "alt", country = "CAN", path = "data/") # Coarse resolution altCAN
## class : RasterLayer ## dimensions : 4992, 10620, 53015040 (nrow, ncol, ncell) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : -141.1, -52.6, 41.6, 83.2 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 ## data source : /Users/mariehbrice/Documents/GitHub/Rspatial/data/CAN_msk_alt.grd ## names : CAN_msk_alt ## values : -108, 5800 (min, max)
# alt90 <- getData('SRTM', lon = -73.7, lat = 45.5) # Fine resolution
As you can see, the resolution of the altitude raster (650x926m) is a lot coarser than that of canopee_mtl (1x1m). We will therefore use the altitude raster for example of raster manipulations to reduce computation time.
Retrieving free GIS data: getData
Retrieve a raster of altitude for Canada.
plot(altCAN)
Exercise
Exercise
use
getData()to retrieve Canadian boundary map at the provincial level?getDataTransform it to a
sfobject usingst_as_sf()Subset the Quebec boundary in the column
NAME_1Reproject your new object using the same projection that of the land use polygons (
st_crs(LU_val)or with the EPSG code32188)Try to plot the geometry of the Quebec polygon.
Try
qc_simple_prj <- st_simplify(qc_prj, dTolerance = 500, preserveTopology = F)and plot the geometry of this new object. Is there a difference?
Exercise
plot(st_geometry(qc_simple_prj))
Geometry manipulations for vector data
Buffer
ruisso_buf <- st_buffer(st_geometry(ruisso_mean), dist = 0.01)
## Warning in st_buffer.sfc(st_geometry(ruisso_mean), dist = 0.01): st_buffer ## does not correctly buffer longitude/latitude data
## dist is assumed to be in decimal degrees (arc_degrees).
dist is assumed to be in decimal degreesThis message indicates that sf assumes a distance value is given in degrees
st_buffer does not correctly buffer longitude/latitude dataThis warning indicates that the result may be incorrect becausest_bufferexpects coordinates in a 2-D, Euclidian space, which is not the case for longitude latitude data. So we should reproject the data onto a projected CRS.
Buffer
plot(st_geometry(ruisso_mean), pch = 19, cex= .8) plot(ruisso_buf, add = T, border= "blue3")
Change projection: st_transform
We can reproject the sample points using a suitable projection that has units of meters. To do this, we will use the projection from our land use polygons.
(ruisso_proj <- st_transform(ruisso_mean, crs = st_crs(LU_val)))
## Simple feature collection with 50 features and 35 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 270615.3 ymin: 5031758 xmax: 304436.2 ymax: 5061940 ## epsg (SRID): 32188 ## proj4string: +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ## # A tibble: 50 x 36 ## sample_pts O2 Conductivite pH Temperature COLI Ag Al ## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 AAO-0.0 8.14 1180. 7.80 17.2 188. 0.100 295. ## 2 AAO-3.3P6 9.50 1378. 7.89 16.4 18599. 0.100 432. ## 3 AAO-3.5 8.79 1281. 7.74 14.3 612. 0.100 150. ## 4 AAO-3.6 9.00 1128. 7.97 16.1 461. 0.100 217. ## 5 AAO-6.4P12 10.1 1461. 7.96 18.2 147. 0.100 33.1 ## 6 AAO-6.5 7.40 567. 8.12 16.5 1219. 0.100 109. ## 7 ADM-1 9.97 333. 8.19 21.2 58.6 0.100 119. ## 8 ANG-2 9.70 399. 8.33 20.2 41.0 0.100 60.7 ## 9 BER-0.0 7.18 803. 7.61 17.1 305. 0.100 140. ## 10 BER-0.7P1 9.25 751. 7.76 17.1 1431. 0.100 75.7 ## # ... with 40 more rows, and 28 more variables: As <dbl>, Ba <dbl>, ## # Be <dbl>, Ca <dbl>, Cd <dbl>, Co <dbl>, COT <dbl>, Cr <dbl>, Cu <dbl>, ## # Fe <dbl>, K <dbl>, Mg <dbl>, Mn <dbl>, Mo <dbl>, Na <dbl>, NH3 <dbl>, ## # Ni <dbl>, Ptot <dbl>, Pb <dbl>, MES <dbl>, Sb <dbl>, Se <dbl>, ## # Sn2 <dbl>, Tl <dbl>, U <dbl>, V <dbl>, Zn <dbl>, geometry <POINT [m]>
Buffer
ruisso_buf <- st_buffer(st_geometry(ruisso_proj), dist = 500) # add an attribute for sample points id ruisso_buf <- st_sf(sample_pts = ruisso_mean$sample_pts, ruisso_buf) par(mar=c(0,0,0,0)) plot(st_geometry(ruisso_proj), pch = 19, cex= .5) plot(st_geometry(ruisso_buf), add = T, border= "blue3")
Exercise
Exercise
Reproject the
courseauobject using the land use projection and name itcourseau_projCreate a 300m buffer around each watercourse and name it
courseau_buf. Plot the resulting geometry withcourseau_projon top.Try
st_union()oncourseau_buf. Plot the resulting geometry and compare withcourseau_buf.Try
st_centroid()onruisso_buf. Plot the resulting geometry withruisso_bufunder.
For more geometric transformation such st_difference(), st_convex_hull(), st_intersection() see sf vignette
Exercise
Exercise
Exercise
Simplify and reclassify land use types
LU_reclas <- LU_val %>%
dplyr::select(ID, UTIL_SOL) %>%
mutate(UTIL_SOL = case_when(UTIL_SOL %in% c(100:114) ~ "resid",
UTIL_SOL %in% c(200:520) ~ "public_build",
UTIL_SOL == 600 ~ "green",
UTIL_SOL %in% c(700:760) ~ "road",
UTIL_SOL == 800 ~ "agri",
UTIL_SOL == 900 ~ "vacant",
UTIL_SOL == 1000 ~ "water",
UTIL_SOL == 1100 ~ "golf"))
Intersecting polygons for area calculations
We now want to find the proportion of area per buffer occupied by different land use types. st_intersection() can be used to obtain the intersection of two geometries (i.e. the area covered by both).
LU_inters <- st_intersection(LU_reclas, ruisso_buf)
## Warning: attribute variables are assumed to be spatially constant ## throughout all geometries
Intersecting polygons for area calculations
plot(LU_inters["UTIL_SOL"])
Intersecting polygons for area calculations
# add in a column and compute area (area of each land use poly in intersect layer) LU_inters$areaLU <- st_area(LU_inters) # group data by sample points and land use type and calculate the total area for each land use per buffer LU_area <- LU_inters %>% group_by(sample_pts, UTIL_SOL) %>% summarise(areaLU = sum(areaLU)) # remove geometry st_geometry(LU_area) <- NULL
Intersecting polygons for area calculations
# buffer with 500m radius so area is pi*500^2 # try st_area(ruisso_buf) # compute proportion for each land use LU_area$propLU <- as.numeric(LU_area$areaLU/(pi*500^2)) # transform to wide format and create a new data frame containing our new landscape variables env_var <- reshape2::dcast(data = LU_area, formula = sample_pts ~ UTIL_SOL, fill = 0)
## Using propLU as value column: use value.var to override.
Geometry manipulations for raster data
Crop a raster for faster computation
crop() will decrease the extent of the raster using the extent of Montreal area.
extMtl <- extent(c(xmin=-73.97342, xmax=-73.47562, ymin=45.39904, ymax=45.73252)) alt_crop <- crop(altCAN, extMtl) # Crop the raster
Crop a raster for faster computation
crop() will decrease the extent of the raster using the extent of Montreal area.
plot(alt_crop)
Reproject a raster
We can use projectRaster() to transform the CRS of one spatial object to match another spatial object
alt_crop
## class : RasterLayer ## dimensions : 40, 60, 2400 (nrow, ncol, ncell) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : -73.975, -73.475, 45.4, 45.73333 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 ## data source : in memory ## names : CAN_msk_alt ## values : 0, 179 (min, max)
alt_proj <- projectRaster(alt_crop, crs = st_crs(LU_val)[[2]]) alt_proj
## class : RasterLayer ## dimensions : 48, 72, 3456 (nrow, ncol, ncell) ## resolution : 650, 926 (x, y) ## extent : 263713, 310513, 5025305, 5069753 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ## data source : in memory ## names : CAN_msk_alt ## values : 0.9513865, 173.2387 (min, max)
Extract and summarise raster values in buffer
Compute the mean altitude per buffer
sample_alt <- raster::extract(alt_proj, as(ruisso_buf, "Spatial"), fun=mean, na.rm=TRUE) head(sample_alt)
## [,1] ## [1,] 23.51085 ## [2,] 28.15839 ## [3,] 28.15839 ## [4,] 28.22916 ## [5,] 32.50139 ## [6,] 30.12374
# Add this new environmental variables to our data frame env_var <- cbind.data.frame(env_var, altitude = sample_alt)
See the package velox for faster raster extraction
Statistical analyses
Now that we have our response variables (mean measurements of water quality in different watercourses) our environmental variables of interest (% of different land use types in a 500m buffer and mean altitude in a 500m buffer), we could performed various statistical analyses. For instance, we could use these data and run a RDA.
Custom map
Custom map
plot() allows combination of multiple layers of information in a same graph, with the argument add = T
plot(canopee_mtl) plot(st_geometry(courseau_proj), col = "#034ebf", lwd = 1.5, add = T) plot(st_geometry(ruisso_proj), pch = 19, col = "#B00C0C", add = T)
Custom map
Change the color palette for the raster levels. There are 5 levels and only 2 are pertinent: 3 = low vegetation cover and 4 = canopy (see here for info on the classification).
myPal <- c("white", "white", "#BAE4B3", "#006D2C", "white")
plot(canopee_mtl, col = myPal, legend = F)
plot(st_geometry(subset(LU_val, UTIL_SOL==1000)), col = "#7eb0fc", border = "#7eb0fc",add=T) # UTIL_SOL==1000 -> rivers
plot(st_geometry(courseau_proj), col = "#034ebf", lwd = 1.5, add = T)
plot(st_geometry(ruisso_proj), pch = 19, col = "#B00C0C", add = T)
legend("bottomright", legend = c("Low vegetation", "Canopy"), fill = c("#BAE4B3", "#006D2C"),
border="white", bty = "n")
Custom map
Add an inset
myPal <- c("white", "white", "#BAE4B3", "#006D2C", "white")
par(fig = c(0,1,0,1))
plot(canopee_mtl, col = myPal, legend = F)
plot(st_geometry(subset(LU_val, UTIL_SOL==1000)), col = "#7eb0fc", border = "#7eb0fc",add=T) # UTIL_SOL==1000 -> rivers
plot(st_geometry(courseau_proj), col = "#034ebf", lwd = 1.5, add = T)
plot(st_geometry(ruisso_proj), pch = 19, col = "#B00C0C", add = T)
legend("bottomright", legend = c("Low vegetation", "Canopy"), fill = c("#BAE4B3", "#006D2C"),
border="white", bty = "n")
par(fig = c(0.08, 0.6, 0.42, 1),xpd=T, new = T)
plot(st_geometry(qc_simple_prj), col = "grey85")
points(286543, 5038936, pch = 17, col= "#B00C0C")
Interactive map with mapview
Interactive map
ruisso_env <- left_join(ruisso_proj, env_var, by = "sample_pts")
## Warning: Column `sample_pts` joining character vector and factor, coercing ## into character vector
Interactive map
mapview(ruisso_env, zcol = 'COLI', legend=TRUE, layer.name = "E. coli density")
Interactive map
mapview(ruisso_env, zcol = "resid", legend = TRUE, layer.name = "% of residential area")
Interactive map
myPal <- colorRampPalette(brewer.pal(9, "YlGnBu")) mapview(ruisso_env, zcol = "resid", cex = "COLI", legend = TRUE, col.regions = myPal, layer.name = "% of residential area")
Ressource for mapping in R
Custom maps
Introduction to visualising spatial data in R
Interactive maps